Dynamics of the Chemical Master Equation, a strip of chains of equations in d-dimensional space. 
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We investigate the multi-chain version of the Chemical Master Equation, when there are transitions between 
different states inside the long chains, as well as transitions between (a few) different chains. In the discrete 
version, such a model can describe the connected diffusion processes with jumps between different types. We 
apply the Hamilton- Jacobi equation to solve some aspects of the model. We derive exact (in the limit of infinite 
number of particles) results for the dynamic of the maximum of the distribution and the variance of distribution. 



PACS numbers: 87.10.-e, 02.50.Ey 



I. INTRODUCTION 



The Chemical Master Equation (CME) 01,11.11,01 is 
the main tool to investigate chemical reactions with few 
molecules. This statistical physics model describes also the 
branching-annihilation processes and (in strip version, con- 
sidered in the current article) the discrete version of the con- 
nected diffusion. 

For reactions with few molecules the number of molecules 
fluctuate, and we can describe the reaction process through 
the probabilities of having a given number of molecules. 
While working with D types of molecules, we work with 
the probabilities of having given integer number of molecules 
P(Xi...Xd, t), thus the variables are located at the nodes of a 
D-dimensional lattice, where integer X{ describe the number 
of molecules of the i-th type. One solves a linear system of 
differential equations for P(Xi...Xd, t) to define the proba- 
bility of having the given numbers of molecules at time t. The 
CME is a linear equation, like the Shroedinger equation for 
the wave function. We assume that the system does not in- 
teract with the environment, and during the measurement we 
observe just some collection of integers for the different types 
of molecules, defined by probabilities P{X\...Xd). 

While the numerical solution of CME is possible using the 
Gilepsee algorithm, the analytical investigation of CME is 
very important. The problem was solved exactly in the infinite 
system size (the maximal number of particles) limit for D = 1 
using quantum mechanics or mapping to the Hamilton-Jacobi 
equation (HJE) |1,|1,|1 Ah. In ^ an im P licit expres- 
sion was obtained for the dynamics of the population distri- 
bution variance, following the HJE method of 01, 01, 01 ■ 
To obtain an exact complete solution of CME in dimensions 
D > 2 is a very hard problem and only a few solutions are 
known, see ITill . 

Sometimes the CME is organized as a strip of systems of 
equations 0-O1- In strips of one dimensional chains 
of equations (the variables are located on 1-D axes and the 
time derivative of some variable depends on its neighbors) are 
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FIG. 1: The CME for the case d = 1,D = 1,K = 1. Possible 
states of the system (indicated as boxes with the number of molecules 
inside) are located on the 1-d space (D — 1). There are transitions 
only between neighbors (K = 1). There is only one chain (d — 1). 
The back transition rates is _R_i and the forward transition rate is 
Ri, see Eq. (1). 



considered as a model of coupled diffusion processes: a situ- 
ation in 2 dimensional space, when there is a diffusion in one 
dimension, and jump process in the other dimension (continu- 
ous diffusion later is replaced by discrete jumps). iTot] consid- 
ers a strip of three chains of equation, each of chains defined 
in 2-dimensional space, identifying the model with genetic 
switch phenomenon. Again we have a diffusion, completed 
by jump process. 

In the current article we formulate the following general 
problem: we investigate the CME, organized as a strip of d 
chains of equations, each defined in D-dimensional space with 
a maximal number of a given type of molecule N, where d -C 
N. 

We will also give differential equations, defining the solu- 
tion of the moments of distributions of the model for the gen- 
eral d, D case. Two mathematical methods have been applied 
to solve the CME in the limit of large number of molecules: 
diffusion approximation and the HJE. The HJE method allows 
to apply the methods of Hamilton mechanics OH, me Hamil- 
ton equation for the characteristics 01, JH- F° r the recent 
introduction to the method see rt2lll . In this article we consider 
the situation when the HJE method is valid while describing 
the solution of the CME. 

Let us first review known results of the CME. In the sim- 
plest case (D = l,d = 1) and one step processes (K = 1), 
we have a system of equations, see Fig. 1 for illustration, 
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FIG. 2: The CME for the case d = 1, D = 1, K = 2. Possible state 
of the system are located on the line (D = 1). We illustrate possible 
transitions from the position n. There are 2-step transitions (K = 2). 
There is only one chain (d — 1). The back 1-step transition rate is 
and the forward 1-step transition rate is R 1 . The 2-step back 
transition rate is R-2 and the two-step forward transition rate is R2 



where < X < N, and N is the maximal possible number 
of molecules, a large number. Ri(X) is the forward 1-step 
transition rate, and R_i(X) is the back 1-step transition rate. 

A similar equations has been formulated for the virus evo- 
lution problems. In case of virus evolution models l2"2tl . l23tl . 
there is inhomogeneous term (fitness function) on the right 
hand side of equation, and P(X) describe the density of dif- 
ferent types of viruses. 

While in virus evolution models one can consider a mixed 
state at the start, in case of the CME, by Eq.(l), there is a 
definite value n for X at the start. Therefore 



P(X,0)=5 X ,n 



(2) 



If during an elementary reaction event there is a consump- 
tion or birth of several molecules with maximal consumption 
(birth) number K, then our system of equations is modified: 



dP(X,t) 
Ndt 



= E 



R m (X -m)P(X -m,t) (3) 



-K<m<K 



Here R± m (X),m ^ are the rates of m-step transi- 
tions, therefore are nonnegative numbers, while Rq(X) = 
— Yl-K<m<K Rm{X) (a generalization of the expression by 
Eq.(l)) ensures a probability conservation condition, see Fig. 
2 for an illustration. The probability rate in the CME are de- 
fined via reaction rates in case of infinite number of molecules, 
see JU. Different rates R-s in our equations describe differ- 
ent elementary reaction events. In the literature, models with 
K < 3 has been considered. The derivative dP(X)/dt de- 
pends on P(X') for the neighbors and that's why we name 
our system of equations as a "chain" of equations. 
Consider the case of D-dimensional space, 



dP(X,t) 
Ndt 



K 



= E E Rr(X - n)P{X - n,t) (4) 



1=1 m = -K 



Here X is a D-dimensional vector, where D is the number of 
different type of molecules participating in the reaction, and 
for ft 7^ the nonnegative number i? r -j describe the probabil- 
ities of the jump with the change of the number of the 1-th 
type via ni, 1 < I < D, and negative number Rq ensures 
the probability conservation. As an illustration consider the 
case D — 4. The process, when the molecules of types 1, 
2, 3 collide giving a birth of 3 molecules of type 4 is de- 
scribed via a coefficient i?-i,-i.-i,3(xi + 1,X2 + 1,^3 + 



1,2:4 — 3). If these three molecules just catalyze the birth of 
one molecule of type 4 the process is described via coefficient 

Ro,0,oAxi,X 2 ,X 3 ,X4 - 1). 

If there are n reactants, then the coefficient in the CME R is 
connected with the coefficient in the chemical kinetic equation 
R with the scaling R = R/N n , see H. 

The Hamilton-Jacobi equation (HJE) is a powerful method 
to investigate the CME at large values of N. For illustration 
consider Eq.(l). For the case Af>lwe assume an ansatz 



P(X,t) = exp[Nu{x,t)\, 
X 

x = — . 

TV 



(5) 



The ansatz by Eq.(5) is similar to the WKB approximation in 
quantum mechanics [9]. Eq.(5) gives with 1/N accuracy: 



P(X±l,t) =exp[Nu(x,t)±u'(x)], 

X 

X= N' 



(6) 



and one obtains |0] : 



du(x, t) 
dt 



+ H(x,u') = 



—H(x,p) = r_(x)e p + r+(x)e p — (r+(x) + r_(x)) 

ox 

where we denote r_(x) = R-\{Nx),r + (x) = R+(Nx). 
Equation (7) is a partial differential equation of a special 
form, the Hamilton-Jacobi equation. According to the stan- 
dard handbook of classical mechanics lfl8ll . it can be solved 
using the solution of Hamilton's equation for x( t), p it), where 
x(t) describes the characteristic curves ofHJETO.lfli.llllll. 
The N in the dominator of the left hand side of Eqs.(3),(4) is 
written arbitrary. We are choosing a proper scaling N n with 
some integer n, to get a smooth, N independent functions r(x) 
at the limit — > 00. We can consider the more general case 
N n , n > 1, then rescaling t — > t/N 11 " 1 return to the case 
of article. The HJE approach allows to define the whole dis- 
tribution (P(X) or u(x)). In the following sections we will 
derive the HJE version of the CME for more involved situa- 
tions than Eq.(3). In this work we are mainly interested in the 
dynamics of the maximum of distribution (the point x where 
u' x (x, t) = 0) and in the variance of the distribution P(X, t), 
given by u"(x, t) at the point of maximum. For our purposes 
there is no need to apply a powerful method of characteristics. 
We will use directly the HJE. We are interested in the solu- 
tion u(x,t) near the maximum of the distribution, x = y(t) 
in our notation. We can derive exact relations by expanding in 
the HJE the solution in degrees of (x — y(t)). Such a method 
has been applied in ifioll . where we derived the Van Kampen's 
large system size expansion method from HJE. 

In Section II-A we derive HJE for the d > 1, D = 1 case. 
While we derive the HJE to describe the full distribution, it 
is more informative to calculate the motion of the maximum 
and the variance. We derive corresponding ODEs in Sections 
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II-B and II-C. In Section II-D we considered an example of 
D = 1, d = 2. In Section III-A we consider the strip of D- 
dimensional equations. In the Appendix we present the equa- 
tions to define the dynamics of the third moment of distribu- 
tion, investigation of evolution insertion-deletion model l24ll . 
In section III-B we give the finite size corrections to the bulk 
solution, following to ifTHl, ll8tl.ll9tl. 



II. THE STRIP OF ONE DIMENSIONAL CHAINS 

A. Derivation of HJE 

Here we consider the CME, when the probability depends 
on two integer coordinates, / and X, thus we work with 
Pi(X, t). There is a probability conservation condition 



Eq. (13) is a consequence of probability conservation con- 
dition during the transitions between different chains at the 
position X, see model by Eq.(31) as an example. 
Let us consider the following ansatz: 



J2 E Pi(x,t) = i 



(8) 



KKd 0<X<N 



In Eq.(8) the maximal number of molecules is constrained by 
N, a reasonable constraint in case of chemical reaction in a 
finite volume, and d <C N. 

In principle only one boundary condition at X = can be 
present or some large parameter N . In both cases HJE method 
can be applied. 

We consider a large N limit: 



N -> oo. 
There is a system of equations 

dPt(X,t) 



(9) 



Ndt 



= £ ]T R lnm (X - m)P n (X - m,t%lO) 



Kn<d —K<m<K 



There are transitions along the states of one chain (maxi- 
mal jump is K-step), described by nonnegative coefficients 
Rum , m ^ and there are jumps between different chains, 
described by nonnegative Ri nm with I ^ n. There are nega- 
tive Ruo, and there is an equation 



J2 J2 Rlnrn(X)=Q 



(ID 



KKd-K<m<K 



to support the probability conservation condition. 

The system is slightly modified near the borders X = 
and X — N. In principle, we can use a one-sided boundary 
condition, but for our case of solution, described by the HJE, 
the boundary conditions are irrelevant. 

Let as assume that at the limit N —> oo (condition a.), the 
rates are described via smooth functions 



X 

X= N 

We consider the models under the condition 

det{A ln } = 0, 



A, 



-K<m<K 



(12) 



(13) 



Pi{X) = v t exp[Nu(x,t)}, 
X 

X= N- 



(14) 



We assume that the maximum of different components 
Pi(X) are near the same X. Thus, if the maximum is at 
the points X\ and X2, then 



\X\ — Xi 
N 



< 1. 



(15) 



The smoothness of functions n„ m in Eq.(12) is crucial for the 
ansatz Eq.(14). In lfl3ll we considered an example, when the 
exponential ansatz p ~ exp[Nu(x)] is invalid for non-smooth 
functions of transition rates. 
We put also the condition b: 



d< N. 

Then we have the following system of equation: 



i-7 



= £ £ ri nm {x)v n e~ 



(16) 



(17) 



Kn<d -K<m<K 



We dropped the terms dvi / dt, assuming d N, where N is 
the system size or some large parameter. 

As the system of equations for vi should be consistent, we 
have the condition: 

det[Mi n (x, u) - qSi n (x)] = 0, 
M ln (x,u')= n nm {x)e- mu ', 

—K<m<K 

<is, 

We assumed that M in {x, u') is a (smooth) analytical func- 
tion of x, and all the terms of M have the same, zero de- 
gree of N (the case of the same, non-zero degree of N could 
be investigated after re-scaling of time). In this case we can 
consider N — > 00 limit in a proper way. We denote x = X/N. 
Eq.(13) gives 



det[M ln (x,0)} =0 



(19) 



Eq. (18) is our main equation. It has a high degree of q = ^ 
and, therefore, it defines a standard HJE with the first degree 
of q as in Eq.(7) and corresponding Hamiltonian implicitly. 
To construct the full solution u(x, t) we have to find such an 
explicit Hamiltonian. Fortunately, there is no need for an ex- 
plicit Hamiltonian to investigate the properties of the solution 
near the maximum of distribution. Expanding the left hand 
side of Eq.(18) in the degrees of q, we re-write Eq.(18) as 



(20) 



Eq.(19) gives: 



H (x,p)\ p=Q = 
du 
^ dx 



(21) 



In investigating the CME we are first interested in the maxi- 
mum of the distribution, then in the variance. 



B. The dynamics of the maximum of distribution 

Let us first consider the dynamics of the maximum. While 
calculating the dynamics of the maximum, it is enough to con- 
sider 0-th and first terms in Eq.(20). 

We assume the ansatz 



u = -^(x-y(t)) 2 . 



(22) 



Let us differentiate Eq.(20) with respect to x, and consider the 
point x — y(t). The higher terms of q l disappear, as q ~ 
(x —y(t)). We obtain: 



dH (x,p) , 
V - | p= o - q x H 1 = 0. 



dp 



(23) 



UsingP = det[Mi n ],H x = —^det[M ln (x,u')-q8 ln ]\ q = , 
we obtain 



dy(t) _ 
dt #i(y,0) [V) ' 



(24) 



Eq.(24) defines the dynamics of concentration in case of the 
infinite number of molecules. It is the first step in the investi- 
gation of the CME. 



C. The dynamics of the variance 

To investigate the dynamics of the variance, we should con- 
sider the q 2 term in Eq.(20). We will consider the solution near 
x = y(t). Dropping the higher terms in (x — y(t)), we obtain: 



dx 



2 [q 2 H 2 -qH x +H )\ P=0 =0. 



(25) 



Using dH ^.' p ^ = H' x + H'p' x , we calculate 



d 2 H(x,p) 
d 2 x 



H' x ' x + 2H' x ' p p' x + H' p p'^ x + H'; p (p' x ) 2 . (26) 



Putting Eq.(26) into Eq.(25), we obtain 

d 2 



dx 2 



[q 2 H 2 - g#i + H ] 



2{q' x ) 2 H 2 - qlJiy - 2q' x [H' lx + H[ p p' x ] 

^OxpPx + H()pPxx + Hoppi 



-Hqxx + ZHoxpP'x + HonPxx + Hq pp (p x ) 2 — 0. (27) 



1,n-1 



2,n-1 





H 

[2^1 



1,11+1 



2,11+1 



FIG. 3: The CME given by Eq.(31) for the case d=2,D = 1,K = 
1. Possible state of the system are located on the two 1-d lines (D = 
1). There are only one-step transitions inside the chains (K — 1). 
There are two chains (d = 2). The backward transition rate is pn 
and the forward transition rate is KN. The transition rate from the 
1-st chain to the second is nh and from the second chain to the first 
is Nf 



Using the ansatz (22), we derive from Eq.(27) 

Vm + V 2 [2b 2 H 2 + 2bH[ p + H'J pp ] + 

+V[-2bH[ x - 2H^ p ] + H? )xx = 0. (28) 

We use the notation b(y) from Eq.(18). Using that Hq xx = 
0, -bH\ - H' 0p = 0, we obtain 

VHx + V 2 [2b 2 H 2 + 2bH[ p + H^ pp ] + 

+V[-2bH[ x - 2H^ xp ] = 0. (29) 

Using dy/dt = b(y), and denoting Q = 1/V, we finally have 

^=A(y) + QB(y) 
dy 



B(y) 



A{y) 



b(y)H 1 (y,0) ' 
[2b 2 H 2 + 2bH' lp + H'J pp ] 
b(y)H x {y,Q) 



(30) 



Eq. (30) defines the dynamic of the variance < (x— y(t)) 2 >- 
Q 

N ■ 



D. The case of two 1-D chains, the genetic switch model 

Consider the case lfl7ll : 
dP u 



dt 

dP 2 ,n 

dt 



fcJV[Pi,„_i - Pi, n ] + p[(n + l)Pi, n+ i - nPi, n ] 

-nhPi, n + NfP 2 . n 
kN[P 2in ^ - P 2 . n ] + p[(n + l)P 3>n+1 - nP 2 . n ] 

+nhP 1 , n -NfP 2 ,m 



Thus, there are back (pn) and forward (KN) transition rates in 
the 1-d chains, as well as nh and Nf transition rates between 
different chains. Fig. 3 illustrates our CME. 



The system is modified at the boundaries: 



dt 
dP 2 ,o 
dt 



= -kNP ha + P P hl +NfP 2fi 
= -kNRo + pRi-fNFzn 



dPi 



N 



dt 
dP 2 ,N 



dt 



= kNP hN -! - P NP hN - NhP^ N + fNP 2 , 
= kNP 2 , N -i - P NP 2tN + NhP^ N - fNP 2t 



N 



N 



(32) 



Eq.(18) gives 

Mn = He 



1) + px(e u - 1) - xh 
M 12 = / 

M 22 = k{e- u ' - 1) + px(e u ' - 1) - / 

M 21 = xh. (33) 

Thus we define the function b(y) for the dynamic of the 
maximum: 

^ = b(y) = k-py (34) 

and functions A, B to define the dynamic of the variance: 

k + py 



My) 



B(y) 



k- py 
2p 



(35) 



k - py 

Our analytical results are in excellent agreement with numer- 
ical solutions of the original system (31) (see Fig. (4,5)). 

Let us derive the time dependence of the ratio at the 
maximum of distribution. From Eq.(18) we have for our case: 



viMu\ x=y ( t ) 
Using Eq. (33) we get: 



v 2 M- 



12\x=y(t) 







''2 , 



\x=y(t) 



y(t)h 



(36) 



(37) 



vi- — f 
The comparison of the numerics with our analytical results is 
given in Fig. 6. While y(t) and Q(t) follow to theoretical re- 
sults immediately, after period of time T <C 1, the v 2 /v\ fol- 
lows to our theoretical result after period of time T <C N. We 
derived y(t), Q(t) using the expressions of P\{X) + P 2 (x), 
thats why y(t),Q(t) follow the theoretical dynamics rather 
faster than v 2 /v\. As y(t), Q(t) do not depend on the ratio 
v 2 (0)/vi (0) of initial values of probabilities at the maximum 
point, we don't give these values in the Figures 4,5. 



III. STRIP OF D-DIMENSIONAL SYSTEMS OF 
EQUATIONS 

A. The bulk solution of the model 

Consider the following system of equations 

dPi(X) 



Ndt 




0.0 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 




1.4 " 

FIG. 4: The model by Eq.(31), N = 100, p = l,h = 0.005, / = 
0.02. The direct solution of the system of equations (smooth line) 
versus the analytical results (solid circles). y(t) describes the point 
of maximum, 1/Q(t) is the variance of distribution. 

y 

0.50 1 




FIG. 5: The model by Eq.(31), N = 4, p = 1, h = 0.005,/ = 
0.02. The direct solution of the system of equations (smooth line) 
versus the analytical results (solid circles). y(t) describes the point 
of maximum, 1/Q(t) is the variance of distribution. 



Here the sum via rh is over the vector space of the stoichiom- 
etry vectors of all of reactions. Now there are transitions in 
the same D-dimensional space with the rate Ru^ and tran- 
sitions between different D-dimensional spaces with the rates 
Rinm, I 7^ n. The negative Ruq are chosen to ensure the prob- 
ability conservation condition. 



53 Rlnm{X) = 



(39) 



KKd 



where at the limit — > oo 



Kn<d rh 



R lnrh (X -rn) = [r lnA (x) + + 0(^)]A fi (40) 




FIG. 6: The model by Eq.(31), N = 100, p = 1, h = 0.005, / = 
0.02. The numerical result for 1)2/ 'vi (smooth lines) for different 
initial values versus the analytical result (dashed line). While the 
numerical results for y(t) and Q(t) actually are identical with the 
corresponding theoretical results, regardless of the initial value of 
i'2(0)/«i(0), the ratio V2(t)/vi(t) coincides with the theoretical 
value after some period of time T, where T <C N. 



where the scaling n is the same for all index sets (7, n, to). We 
can rescale the time, then putting n = 0. 
Assuming an ansatz 



Pi(X,t) = viexp[Nu(2,t)], 
we have the following system of equations: 



Vi 



du 
dt 



E 



Kn<d rh 



■ E fc rn k u' k 



(41) 



(42) 



For the consistency of the system we again assume that 

det[Mi n (x, u') - qdin] = 0, 



du 



(43) 



We assume an ansatz 



u(x, t) 



2 ^ 

id 



Vy(au-W(*))0ci -!&(«)) (44) 



As before, we consider the decomposition of the determinant 
via Hi,0 < I < d. To calculate the dynamic of the maxi- 
mum, we neglect the (—q) terms for l-s higher than 1. Let us 
differentiate the expression Hq — qH\ with respect to Xk- 



Differentiating the equation Hq — qH\ + H2q 2 with respect 
to Xi,Xj , we obtain for the correlation functions: 

dVij[t) 



l n 



+ E V ^ E H ^n Vjn + E V ifo E H 'lPn V in 
In In 

ill 

1 1 

E H o*w v » E u '** m v * = °- ^ 
1 1 

This is a complete system of equations. 

B. The finite N corrections 

Let as assume that the finite size corrections are u\/N for 
u and hi/N for vi. Thus, we have: 



Pl(X,t) = {vi + -±)exp[Nu(x,t) +u 1 {x,t)} (48) 

The Hamilton-Jacobi equation for the correction terms will 
be: 

dui , du 



QinVn + ^n n , n {x)h n e £* m * w i 

1 < n < d. n,m 

Qin = y^(r lniA (x) + r lnrn (x) * G) 

111 

1 D 



2 dxidx k 

i,k—l 



d u ^ dui 

rrnmk - } m k - 



dxi 



(49) 



From the consistency of vi system we get the following 
equation for u[ : 



du 

det[Qln ftffyn] = 



(50) 



Having the solution for u[ , we define the correction terms h r , 



E 



dH a (x,p) dpi 
dpi dx k 



l' Xk Hi = 0. 



(45) 



^ = -Vki and q' Xk = £. V ki ^, we have 



dpi 

dyjjt) 
dt 



dt 

dH (x,p) I 
dpi \P=° 

Hi(y,0) ' 



(46) 



C. When we need in several Hi ? 

In l28ll has been investigated the D = 2, d = 2 model in 
our classification. 



dP„ 



dt 



{A - f(n))P mn + g(n)Q 



dQ 

III 

dt 



f(n)P mn + {A + a{E^ - 1) - g(n))Q mn (51) 
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where E{f(n) = f(n + j) and A = (££ - l)n + 7^ - 
l)m + 7&?7i(£r 1 - 1). 

The authors applied HJE approach to investigate the steady 
state, using same Hamiltonian, proportional to Hq in our noti- 
fication. Their Hamiltonian gives exact equation for the steady 
state, while for the dynamic it is only an approximation, and 
to get correct dynamics in l28ll has been done some rescaling. 

Our approach provides the exact dynamics in this case, us- 
ing Hi for the dynamics of the maximum, and H2 for the 
dynamics of the variance. 

For Hi(p = 0) we obtain an expression: 

H 1 (p = 0,n/N) = f(n)+g(n) (52) 

It is a nontrivial function of n/N, therefore Ho alone is not 
enough to give a correct dynamics of the maximum. 

IV. CONCLUSION 

In the present article we solved the Chemical Master Equa- 
tion for the case of a strip of 1 dimensional chains of equa- 
tions. Such a master equation appears in the problem of con- 
nected diffusion in discrete time: diffusion in space direction 
plus random jumping between discrete set of states. We found 
a solution for the general case of the model when the width of 
the strip is smaller than the length, and calculated the dynam- 
ics of the maximum and the variance in the large size limit. 
We applied our method to get the solution for the problem of a 
gene switch model. The numerics confirmers that the dynam- 
ics of the maximum and the variance follows to our analytical 
results after very short period of time. 

We also investigated the case of a d-strip of equation chains, 
each one defined in D-dimensional space, and gave a sys- 
tem of ODEs to define the motion of the maximum point and 
the dynamics of the variance. In the appendix we give the 
ODE systems for the investigation of the insertion-deletion- 
base substitution problem in evolution. 

We considered the situation when Hamilton Jacobi equa- 
tion is valid, assuming that the maximum of distribution for 
different chains are at the same point. This was the case in 



the example considered in section II-D, satisfying two con- 
ditions: a) where the rates r(x) are smooth functions and b) 
d <C N. In [24] a D = 2 model has been considered where 
condition b) is broken. As a result, the solution of the sys- 
tem of equations could not be described via the HJE ansatz 
of a smooth function u(x\, X2). To investigate the validity of 
our main assumption Eq.(14), one needs to undertake further 
work with more thorough mathematics. Because nowadays 
the HJE is becoming popular in solving the CME, as well as 
for the solution of the diffusion process in case of weak noise, 
this issue (the validity of the HJE ansatz) is becoming more 
important. Another related fundamental open problem is the 
limit of the application of diffusion approximation in master 
equations of evolution research lf25ll . While the diffusion ap- 
proach is a rather pure ap prox imation in case of the inhomo- 
geneous master equation ll ll Il3ll . it is a very closely related 
technique to solve the HJE for homogenous master equations 
(as considered in the current article), see also rt2lll . Clarifying 
the limits of the powerful HJE method could help to clarify 
the limits of the diffusion method as well. 

Besides the HJE (WKB in literature) method, quantum me- 
chanical approach and the moment closing method have been 
applied to investigate CME. The advantage of the method: 
it allows a series expansion in higher degrees of 1/N, a uni- 
form convergence of our approximation, which lucks for mo- 
ment closing approach. Sometimes a quantum mechanical ap- 
proach is as effective as the HJE one, but HJE has wider area 
of applications. 

Let us briefly discuss some related works. In ll26ll has been 
solve the related model, using the generating function method, 
and has been derived several impressive results in the large 
number of particles limit. This exact method is heavier than 
HJE approach, used in current work, and has narrow area of 
application. In f27ll the gene expression has been investigated 
using spectral method. Perhaps the method is good for the 
numerics. 
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Appendix B: High dimensional case 

1. The dynamics in D-dimensional space 

For the model given by Eq.(4) we get the HJE for the func- 
tion u(x) with following Hamiltonian: 

du . 

— + H(x,f) = 0, 

-H(S,0) = [Rft(Nx-n)e-^-R H (Nx)] . (B.l) 

ft 

Let us consider an ansatz 



Appendix A: Case of two 1-d chains, section II-D 

We consider the model given by Eq. (31). We obtain 
the following analytical solutions for equations Eq.(17) and 
Eq.(24): 



Substituting into the above equation, and differentiating by Xi, 
we get: 



3 3 



The last system of equation is consistent at: 



(B.3) 



y{t) = - + {y --)e-?\ 
P P 



(A.l) 



d Vj _ rr, 



(B.4) 



As a result, we have the following system for the Vij 



Q(y) = y-yo{ 



yo - k/p> 



(A.2) 



where yo is the initial point of the maximum of distribution. 

It is also possible to calculate the variance with more preci- 
sion if we find higher coefficients for the ansatz of u(x,t). To 
find T(y) we use the following equation: 



d 3 [q 2 H 2 -qH 1 +H } 
dx 3 



(A3) 



Using the expressions for the higher derivatives of Hq and Hi , 
we eventually get the following differential equation for T(y) 

at x = y(t): 



dVij Q2 H 

dt dpidpm 1 ™ J 

Ira 

£f dpidx t 1 Z-j dpidxj 



(B.5) 



2. The dynamics for the evolution model with insertions and 
deletions 



Consider the following model, describing the insertions and 
deletions lfl3ll . P{N,L) is the density of the states having 
genome length N and L of + spins, a is the mutation rate, 
b is the deletion rate, c is the rate to insert + or - spins: 



dT 
dy 



bH, 



-T[6Vb 2 - 3bH' lx + 6VbH' lp - 3H^ p + WH^] 
-^b(H[ x - VH[ p 2Vb) 



+Wb{V 2 H'; pp ■ 
-3VH 0xxp — 3V H 0xpp 



V 3 H' " ppp (AA) 



where V = 1/Q, b = dy{t)/dt. 



dp(N,L) = 
dt 

( N N + 1 

+ a(±p(N,L-l) + p(N,L + l)^ 



+ b \p(N + 1, L + 1) +p(N + 1, L)\ 



N + 1 
N 



(B.6) 



C -(p(N-l,L-l)± + p { N-l,L)^ 
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We should re-write the equation for P(N, L) = p(N, L) ( 



Thus, we get the following expression: 



dP(N, L) 
dt 



/ N N +1 

P { N,L)( K N -(a + b )--c^- 

a { ' N -„ L+1 P(N, L—l) + ^P(N, L + l) 



N 



N 



[ ^±lp(N+l,L + l) + - A L + 1 P(N+l,L) 



Nn 



Nn 



*(w p ( N - 1 > L -v + w p t N - 1 >Q 



We have HJE 

- H(x 1 ,x 2 ,pi,p 2 ) =N -(a+b+ c)x\ + 
+a[(x 1 -x 2 )e~ P2 + x 2 e P2 } 
+b[x 2 e Pl+P2 + {x x - x 2 )e Pl ] 

+ ^[x ie - pi - p2 + x ie - pi }. 



(B.7) 



(B.8) 



We are again interested to find the dynamics of the maximum 
and the variance. Let us introduce the following ansatz: 

u(xi,x 2 ,t) = --Vn(xi - yi(t)) 2 - 
-\v 22 {x 2 - y 2 (t)f - V 12 ( Xl - yi(t))(x 2 - y 2 (t)).(B.9) 
For the maximum dynamics we get: 



V2(t) = y2£ e -(2a+c)t + 1[1 

m(t) yw 2 

lim 



s -(2o+c)ij 
V2(t) 1 



mil ^r- 

t^oo yi (t) 



(B.12) 



For the variances we eventually get the following system of 
equations: 



li 



dt 



a[2V 12 + yi V? 2 ] 



+b[y 2 (Vn + V 12 ) 2 - 2V U + ( yi - y 2 )V y 



li J 



':[4Vn + 2V 12 + Vl (V u + V 12 f + yi Vft]. (B.13) 



dt 



a[V 22 ~ 2V 12 + Vl V 12 V 22 } 



+b[-2V 12 + y 2 (V 12 + V 22 )(V 12 + V n ) 

+{y\ - y 2 )VnVi 2 ] 
-|[2^i2 + V 22 + yi(V 12 + V 22 )(V 12 + V n ) 



+yxV n V 12 }. (B.14) 



-a= [c - b)yi 

^-=a(y 1 -2 y2 )-by 2 + ^y 1 . (B.10) 
The solutions of the latter equation is: 

yi(t)=y 10 e^ t 

(B.ll) 



dV; 



2 =a{-W 22 + yi V 2 2 2 ] 



dt 

+b[-2V 22 + y 2 {V 12 + V 22 ) 2 + { m - m )V 2 2 } 

+ ^[yi(V 12 + V 22 ) 2 + yi V 1 2 2 }. (B.15) 



